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1 Introduction 



The Hamiltonian of a system of two scalar particles with mass to in a local potential is in 
the relativistic situation given by 



We are interested in solving the Schrodinger equation for binding and confining potentials, 
that are associated with particle exchange. 

There are several methods for obtaining bound states. Probably the simplest way is to 
directly solve the differential equation in configuration space. Unfortunately we cannot use 
this method in the relativistic case, since we have to deal with a term that contains a square 
root of the differential operator. 

Another possible method is to directly discretize the integral equation in momentum 
space. This comes down to solving an equation of the form 



A third method makes use of the variational principle for the energy and expands the eigen- 
states in terms of basis functions (tp — &i<t>i)- It turns out to be advantageous to be able 
to calculate the kinetic part of the Hamiltonian in momentum space and the potential part 
in configuration space. 

In this report we will calculate the bound states with both methods. The issues we 
concentrate on are ease of numerical implementation, accuracy and stability. Both methods 
are less easy to implement and calculations take up more time compared to the method of 
directly solving the differential equation. To compare the accuracy of the two methods, we 
first need to know whether our codes are correct. To check our codes we have made use of 
several potentials for which the bound state energies and wave functions are known explicitly 
in the nonrelativistic case. In section 2 we will present these potentials. In the next section 
we will look in more detail at the numerical implementation of the methods. In section 4 
we will give some numerical results in the nonrelativistic case and in the next section we 
will look at some physical interesting cases. Finally we will give some conclusions in the last 
section. 



H = 2^/m 2 +p 2 + V(r) 



(1) 




(2) 



2 Exact solutions of potentials in the nonrelativistic case 



There are several potentials for which the bound state energies and wave functions are known 
in the nonrelativistic situation. We have used three of these potentials to check our codes. 
They are the Coulomb potential, the linear potential and the Hulthen potential. Note that 
two of these potentials have some features in common with the Yukawa potential: the 
Coulomb and Hulthen potentials behave like 1/r at r — > 0; the Hulthen potential falls of 
like an exponential as r — > oo, almost the same as the Yukawa potential 

V(r) = ^e-^ r . (3) 

By looking at this potential, Eq. ([|), one realizes immediately that putting /iy = reduces 
the Yukawa potential to the Coulomb potential. The energies of the bound states in the 
nonrelativistic Coulomb case are known exactly and given by 

w - mpI d\ 

2.1 The linear potential 

In the case of the linear potential (V(r) — err) the eigenvalues are exactly known for I = in 
the nonrelativistic case and can be expressed in terms of the zeros of the Airy function Qj. 
The eigenvalues are given by 

ct2X 1/3 



E n = \- X n (5) 

where a is the strength of the potential and A„ is the n-th zero of the Airy function (see 
table |). 

2.2 The Hulthen potential 

The Hulthen potential has the following form 



e 



V(r)=p- , (6) 

where p is the strength of the potential. The wave functions are known for this potential in 
the nonrelativistic case (I = 0) and are given by 

rip(r) = u k (r) = e- a " r {\ - e -^ p ^lmA)^ _ 2e -^ r ), (7) 

where is defined as 

ak = 2(fe+l)MH ' (8) 
p^2Q fc /^H,i) j g a j aco ^,j polynomial ?? and the bound state energies are given by 

E k = (9) 
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3 Implementation 



There are different methods to calculate the bound state energies and wave functions of a 
system. The methods we present here are chosen for different reasons; ease of implementation, 
accuracy, speed, stability etc. First we will present the different approaches to calculate bound 
states. Then we will look in more detail at the methods we use and pay attention to the 
stability and accuracy of these methods. 



3.1 Different approaches to calculate bound states 

The first method we mention is the method in which the differential equation is solved step 
by step. It is very easy to implement. The calculations do not take much time, but it is not 
very accurate. And unfortunately this method gives problems in the relativistic case. 

Discretizing the integral equation is a more accurate method. It is less easy to implement 
and calculations done using this method take up much more time. This method is not 
very well suited to all potentials, since all calculations are done in the momentum space 
representation and we may have to deal with singularities in the potential matrix elements. 

The third method we present expands the eigenfunctions in terms of basis functions and 
does not have this problem. The two parts of the Hamiltonian, the kinetic and potential 
energies respectively, are calculated in either momentum space or configuration space, what- 
ever is more convenient. This method also gives a better accuracy than the first method, but 
it also takes more time to do the calculations. 



3.1.1 Directly solving the differential equation 

Probably the simplest way to determine the bound states of a nonrelativistic system is to 
directly solve the differential equation. We have done this by reducing the Schrodinger 
equation to a set of coupled first order differential equations and solve them using the Runge- 
Kutta method [|l]. This method is only simple to implement for a nonrelativistic system. 
For a relativistic system it becomes much more difficult or even impossible to calculate the 
bound states, since we have to deal with a term that contains a square root of the differential 
operator. 



3.1.2 Discretized integral equation 

The radial part of the Schrodinger equation in spherical coordinates is given by 
h 2 d 2 h 2 l(l + 1) „. A , . . 

— -T2+ 2 L +V(r))Mr) = EMr), (10) 
m ar mr J 

where I is the orbital angular momentum quantum number and ipi(r) is the radial part of 
the total wave function which is defined as 

*(f) = Mr)Y lm {f) = ^-Y lm (r). (11) 

We can also write Eq. ( |To| ) in the momentum space representation. We then find the integral 
equation 

2 /*oo 

Z-fa(p)+ dqq 2 Vi(p,q)^i(q) = Efa(p). (12) 
m Jo 

Here ipi (p) is the Fourier transform of the function ipi (r) and Vi (p, q) is the Fourier trans- 
formed potential. 
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Since this integral equation cannot be solved in closed form, we apply a suitable Gaussian 
integration rule, with abscissas qj and weigths Wj and approximate the integral by a finite 
sum. We first transform the integration variable q to a variable x, so that the integration 
interval changes from [0, oo) to [—1, 1). This is done by the transformation 

q{x) = c\±^, (13) 
1 — x 

where C is a positive scaling parameter. Because of this transformation, the integrand 
contains an additional factor, the Jacobian dq(x)/dx. 

Now that we have approximated the integral by this finite sum, we can take the momentum 
variable p equal to the abscissas. By doing this we obtain an eigenvalue equation in matrix 
form, which has the following form when we choose a total of N abscissas: 

N 

J2 K vM<lj) = Ei>i{qi), (14) 
j=i 



where the matrix-elements of the matrix K are given by: 

q 2 

K i:j = -LSij + W 3 q 2 Vi{q z , qj) (15) 



with Wj = (dq(xj)/dxj)wj. 

The matrix equation ( |14|) would be easier to solve if the matrix would be symmetric. Left- 
multiplying the matrix equation with a diagonal matrix D, with elements Dij = qjSij j \Jw~j, 
and inserting the identity D~ 1 D between the matrix K and the column vector ipi gives us a 
symmetric matrix K = DKD^ 1 . The matrix elements are given by 

q 2 



ni 

for the nonrclativistic situation, and 



Kij = —S i: j + y/wiWjVi(qi, qj) (16) 



Kij = 2 J q 2 + m 2 Sij + y/wiWjVi(qi,qj) (17) 

for the relativistic situation, with vi(p, q) — pqVi(p, q). The matrix equation we solve in this 
approach is given by 

N 

Y,K ij u l {q j ) = Eu l {q i ), (18) 

where ui is the eigenvector, whose j-th element is given by yjw~jqjvpi{qj). 

3.1.3 Expanding the eigenfunctions in terms of basis functions 

In the third approach the eigenvalues arc calculated through expanding the eigenfunctions 
in terms of basis functions. The basis we use in this approach is the one first introduced by 
Olsson ]3j and Weniger Q . The basis functions are known in this basis in both configuration 
and momentum space. In configuration space they are given by 

<l>k,l (r) = N H (2 M r)< e-^Lf +2) (2 M r) , (19) 
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where the normalization is given by 



_ fc!(2 M )3 

y (fc + + 2)! (20) 

and Lf +2) (2H 

is an associated Laguerre polynomial. In the momentum space representa- 
tion the basis functions are given by 

7 I \ at i f V Y +2 D ('+3/2,i+i/2) (p 2 - M 2 \ 
with P k a '^\x) a Jacobi polynomial and the normalization constant is 



- 2y M fc!(fc + 2/ + 2)! 

^ - r(fc + / + 3/2) ■ (22) 

The first five wave functions are plotted in Fig. ^in both configuration and momentum space. 
Note that these functions differ from the familiar Coulomb and harmonic oscillator basis func- 
tions. Ref. [@ gives the mathematical details including the proof of completeness of the basis. 

With these basis functions the matrix element of the potential between the initial and 
final states is calculated in configuration space. As an example we will calculate the matrix 
elements for the Yukawa potential. Other potentials are treated in a similar way. The angular 
part of this matrix element reduces to the product of two Kronecker delta functions. The 
more interesting part for us is the part that depends on the radial coordinates; which looks 
like 

dr<Jr)^h^) UMj(r) . (23 ) 

The integrand can be rewritten in a somewhat different form. This is done by first making 
the following substitution for the reduced radial wave function 

u M {r) = N kl L[ 2l+2) (2^y +1 e-^ = 6 H (r)e-^, (24) 

where u k i(r) is a polynomial. Note that the (reduced) radial wave functions behave like r l+1 . 
When we make this substitution, we find for Eq. ( p3[) 

» 1 

dr-exp(-(2^ + MY)r) u* kil .(r) u kjh {r). (25) 
An integral of the form 

dx e-*f(x), (26) 

where f(x) is a polynomial, can be estimated to the order of machine precision using Gauss- 
Laguerre integration. The expression in Eq. ( |25| ) can be written in the desired form by first 
making the substitution x = (2/i + /iy)r. Discretizing the equation and then making another 
substitution ta — Xj/(2/i + fiy) gives the desired form: 

^w k —u* kih {r k ) u kjh (r k ), (27) 



5 



where Wk is given by 

w k = 7T— • ( 28 ) 

2(1 + fly 

The abscissas Xk and weights Wk are obtained according to the Gauss-Laguerre rule. 

To calculate the bound states of the Hamiltonian given by Eq. ([!]) (or its nonrelativistic 
equivalent), we expand the eigenfunctions in terms of the basis functions, ip = a^i, of 
the Olsson basis and rewrite the eigenvalue equation as 

"</".; = J2 ^Wj) °-i = ECL *- ( 29 ) 



3.2 Determining the accuracy of our methods 

There are several factors that influence the accuracy of our results. First we should make 
a distinction between the accuracy obtained through making a mathematical approximation 
and the accuracy obtained because of numerical reasons. The integrals we need when calcu- 
lating the matrix elements of the Yukawa potential can be calculated exactly. However, all 
other integrals are estimated with a mathematical approximation. 

Of course the accuracy of our results depends on more than only the accuracy with which 
we can calculate the integrals. We have looked at the factors that have an influence on the 
accuracy when we use our second method (discretizing the integral equation). Here we have a 
'scaling' parameter C, which has an influence on the accuracy. Another factor that influences 
the accuracy is the number of abscissas we take in our calculation. 

In table [l] we give the calculated bound states for one case: m = 1.0, p = —2.0, (iy = 0.02 
and a scaling parameter C = 1.1. We have calculated the bound states for five different 
numbers of abscissas. From table [l] we see that the bound states converge. 



Table 1: Bound state energies for the Yukawa potential with m = 1.0, p — —2.0 and /iy = 0.02 
calculated with different numbers of abscissas. 



k 


20 absc. 


50 absc. 


80 absc. 


125 absc. 


155 absc. 





-1.189770 


-0.998801 


-0.973828 


-0.964817 


-0.962900 


1 


-0.279598 


-0.221291 


-0.215165 


-0.213154 


-0.212748 


2 


-0.111285 


-0.079688 


-0.077067 


-0.076322 


-0.076185 


3 


-0.053701 


-0.032550 


-0.031212 


-0.030875 


-0.030818 


4 


-0.028074 


-0.012991 


-0.012274 


-0.012113 


-0.012086 


5 


-0.020798 


-0.004369 


-0.004007 


-0.003934 


-0.003923 



To see what kind of influence the scaling parameter C has on the bound states, we have 
calculated the bound states at different C for the case m — 1.0, p = —2.0 and fiy = 0.02 and 
with 155 abcissas. When we look at the results we see that there is a small domain of scaling 
parameters where the bound states are stable. For scaling parameters that lie outside this 
stable domain, the calculated bound state energies are larger. In all our calculations using 
our second method (discretizing the integral equation) we have taken a scaling parameter of 
C = 1.1, which lies in the middle of the stable domain we found. 

In a similar way we have looked at the factors that have an influence on the accuracy 
when we use our third method (expanding the eigenfunctions in terms of basis functions). 
For this method the number of basis functions we have taken into account has an influence 
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on the accuracy. The parameter p is completely free and can be used to tune the basis 
to the Hamiltonian at hand. Although this 'scaling' parameter can thus be considered as 
a variational one, in practice one should be careful if the integrals one needs to calculate 
are obtained numerically. Then it turns out that only in a limited range of values p is of 
genuinely variational character. 

There are several potentials that we have used and for all these potentials we have looked 
what kind of influence the variables mentioned above have on the results. In this report we 
will only show it for the Yukawa potential. 

We have looked at the influence of the number of basis functions we take into account. 
This was done for the following choice of parameters, which are given below: 

m = 1.0, p = -2.0, p Y = 0.02, p = 0.35 

We have made an expansion of the eigenfunctions in terms of the first 5, 10, 15 or 20 basis 
functions; these eigenvalues are given in table |^. We have also calculated the wave functions 
in both configuration and momentum space, these are given in resp. Figs. | and §. 

Table 2: Bound state energies for the Yukawa potential with m = 1.0, p — —2.0 and p\ = 0.02 
calculated with different numbers of basis functions. 



k 


5 basis func. 


10 basis func. 


15 basis func. 


20 basis func. 





-0.9419469833 


-0.9605381489 


-0.9605921507 


-0.9605921507 


1 


-0.2122372389 


-0.2122966051 


-0.2122966051 


-0.2122966051 


2 


-0.07603943348 


-0.07604002953 


-0.07604002953 


-0.07604002953 


3 


-0.02350080013 


-0.03075730801 


-0.03075850010 


-0.03075850010 


4 




-0.01097238064 


-0.01205313206 


-0.01206004620 


5 






-0.002824902534 


-0.003842115402 



Looking at the values of the bound state energies given in table g, we see that these 
converge. Figs. ^ and || show that the wave functions are smoother when we take more basis 
functions into account. From this evidence we conclude that our results can not be taken 
to be reliable when only a small number of basis functions is used in making the expansion. 
When comparing the results of the expansion to the other methods we used 20 basis functions. 

The 'scaling parameter' p appears in the basis, but it is not a real variational parameter 
since we cannot calculate the matrix elements of the kinetic energy without making math- 
ematical approximations. For several values of py in the range between and 1 we have 
looked at the behaviour of p. When we look at the values of the ground state energy at one 
value of py and different /i, we find the same value (for at least the first six decimal places) 
of this ground state for p > 0.30. For p < 0.30, we find that the ground state is shifted 
upwards and has the 'wrong' value. We also have determined a stable domain in p where the 
first five eigenvalues are the same up to at least six decimal places. For all values of Py the 
stable domain starts at p — 0.30 and ends at some value say /Lt max . These maximum values 
are given in table for several values of /iy ■ 



Table 3: Maximum value of p for a stable domain at different values of py- 



MY 


0.0 


0.001 


0.002 


0.005 


0.01 


0.02 


0.05 


0.1 


0.2 


0.5 


1.0 


Mmax 


0.414 


0.545 


0.520 


0.542 


0.476 


0.367 


0.467 


0.663 


1.0 


1.0 


1.0 



7 



In most calculations we have done, we have taken a value of /i = 0.35 for the 'scaling 
parameter', which lies in the stable domain we have determined above. 

All together we see that both methods (discretizing the integral equation and making an 
expansion in basis functions) are not as easy to implement as the method of directly solving 
the differential equation, although it is not extremely difficult either. Both methods have a 
stable domain in which we find accurate solutions. These domains are different for different 
potentials. The accuracy of both methods is better than the accuracy reached with directly 
solving the differential equation. Comparing the two methods, we see that the results found 
by making an expansion in basis functions are more accurate. 

4 Check of the numerical results in the nonrelativistic 

case 

We have first calculated the bound state energies and wave functions in the nonrelativistic 
situation for some potentials for which we know the solutions. After this we have calculated 
the bound states of the Yukawa potential in the nonrelativistic case. 

4.1 The Coulomb potential 

We have calculated the eigenvalues for the Coulomb potential only with the third method. 
The reason is that the long-range character of the Coulomb potential leads to a singularity 
in momentum space that is difficult to handle numerically and in configuration space to an 
uncertainty whether the asymptotic region is reached in the differential equation. 

We limited ourselves to making an expansion of the eigenfunctions in 20 basis states and 
checked that the first five eigenvalues are accurate. In the calculation we have taken the mass 
m = 1.0, the strength of the potential p = —2.0 and the 'scaling parameter' fi — 0.35. 

The calculated eigenvalues and the exact eigenvalues of the Coulomb potential are given 
in table ^ for the first eight eigenvalues. 

Table 4: The calculated eigenvalues (I = 0) of the Yukawa potential with /iy = and the 
exact eigenvalues of the Coulomb potential. 



n 


exact 


calculated 


1 


-1.0000000000 


-1.0000000000 


2 


-0.2500000000 


-0.2500000000 


3 


-0.1111111111 


-0.1111111641 


4 


-0.0625000000 


-0.0625000000 


5 


-0.0400000000 


-0.0399999619 


6 


-0.0277777777 


-0.0277774334 


7 


-0.0204081630 


-0.0202956199 


8 


-0.0156250000 


-0.0137485265 



4.2 The linear potential 

We have calculated the bound states for the linear potential only with the method of ex- 
panding the eigenfunctions in terms of the Olsson basis functions. The bound states are 
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calculated in the nonrelativistic situation for m = 1.0, a = 1.0, p = 0.90900 and are given in 
table |. 

Since a = 1.0 and m — 1.0 the bound states should be equal to the zeros of the Airy 
function. Comparing the results in table || with the zeros of the Airy function, we see that the 
first excited state is accurate to the seventh decimal. Higher excited states are less accurate; 
the accuracy decreases to an accuracy of ~ 1 : 8000 for the fourth excited state. 

Table 5: Zeros of the Airy function and eigenvalues for the nonrelativistic Schrodinger equa- 
tion with m = 1, a = 1. The parameter p in the basis has the value p = 0.909. 



n 


Airy zero 


eigenvalue 


1 


2.33810741 


2.33810741 


2 


4.08794944 


4.08794976 


3 


5.52055983 


5.52057237 


4 


6.78670809 


6.78687212 


5 


7.94413359 


7.94534443 



In this case we have also calculated the wave functions. These are plotted in Fig. ^. 
4.3 The Hulthen potential 

For the Hulthen potential we have calculated the bound state energies with our third method. 
In table H the calculated and exact bound states are given for the case m = 1.0, p = — 0.5 and 
/iH = 0.1. As 'scaling' parameter we have used here p = 1.0 and we have made an expansion 
in terms of the first 20 basis states. 

Table 6: Calculated and exact bound state energies for the Hulthen potential for m = 1.0, 
p = —0.5 and pn = 0.1. 



k 


calculated 


exact 





-6.002500057 


-6.002500000 


1 


-1.322499990 


-1.322500000 


2 


-0.466944456 


-0.466944444 


3 


-0.180624962 


-0.180625000 


4 


-0.062495589 


-0.062500000 


5 


-0.010959864 


-0.013611111 



We see that the first five calculated and exact bound states agree to at least four decimal 
places. Only the very weakly bound fifth excited state cannot be determined accurately. 

In Fig. [5] the exact wave functions and the wave functions calculated by making an 
expansion into 20 basis states are plotted in configuration space. We see that the wave 
functions we have calculated are the same as the exact wave functions. In Fig. [|the calculated 
wave functions are plotted in momentum space. 

We have also calculated the wave functions by directly solving the differential equation 
(Fig. 0). Here we have used the bound state energies calculated with the third approach as 
an input. 

Note that the wave functions in Fig. fj] arc not normalized. To see whether these wave 
functions are the same as the exact and calculated ones in Fig. [|we can look at the location 
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of the nodes and of the extrema. In table |7] the nodes and extrema are given for the exact 
wave functions and the wave functions calculated by directly solving the differential equation. 
Looking at table (?] we see that the nodes and the maxima coincide. 

Table 7: Nodes and extrema of the exact wave functions and the calculated wave functions 
by directly solving the differential equation for the Hulthen potential. 





nodes 


extrema 


k 


exact 


calculated 


exact 


calculated 









0.400 


0.400 


1 


0.8004 


0.800 - 0.804 


0.304 


0.304 








2.098 


2.098 


2 


0.7616 


0.760 - 0.764 


0.296 


0.296 




2.8514 


2.848 - 2.852 


1.680 


1.680 








5.278 


5.278 


3 


0.7501 


0.750 - 0.752 


0.292 


0.292 




2.6632 


2.660 - 2.664 


1.608 


1.608 




6.3132 


6.312 - 6.316 


4.314 


4.312 - 4.316 








10.146 


10.140 - 10.148 



4.4 The Yukawa potential 

In the case of the Yukawa potential (with py ^ in Eq. (g)) there are no exact solutions 
known. We do know, however, that the eigenvalues in this case are greater than the corre- 
sponding eigenvalues in the Coulomb case. We also know that the eigenvalues should vary 
smoothy when py is varied. 

To find the stable domain in which we want to calculate the eigenstates, we have varied p 
for different values of py. In all calculations we have used a value of p = 0.35, which is in the 
stable domain for all values of py. Now that we know where the stable domain lies, we have 
calculated the eigenstates for the Yukawa potential for several values of py . In every case we 
have used 20 basis states in the expansion and expect the first five states to be accurate. All 
the eigenvalues are calculated for a mass to = 1.0 and a strength of the potential p = —2.0. 
The results are given in table ||. 

We see that for small py there are several bound states. As py increases, the number of 
bound states decreases. For p = —2.0 and py > 1.8 there are no bound states found at all. 
The calculated spectra are also plotted in Fig. ||. In these plots we see very clearly that all 
energies shift upwards, when py increases. 

We have not only looked at the eigenvalues, but also at the wave functions. For one case 
we give here the eigenvalues and wave functions as we have calculated them with the different 
methods. The calculated eigenvalues for the case 

m = 1.0, p= -1.0, py = 0.02 (30) 

are given in table §. Note that only the bound states calculated with the second and third 
method are given in this table. For the first method (directly solving the differential equation) 
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Table 8: Eigenvalues (I — 0) of the Yukawa potential for different values of fiy in the 
nonrelativistic case m = 1.0 and p = —2.0. (a ' ' means a bound state has not been found) 



k 


My = 0.01 


p Y = 0.02 


p Y = 0.05 


/i Y = 0.1 





-0.9801490307 


-0.9605921507 


-0.9036328793 


-0.8141160011 


1 


-0.2305865288 


-0.2122966051 


-0.1635423899 


-0.0998564959 


2 


-0.0923976898 


-0.0760400295 


-0.0387051106 


-0.0064160824 


3 


-0.0447121859 


-0.0307585001 


-0.0061832666 




4 


-0.0233221054 


-0.0120600462 







k 


/i Y = 0.2 


/iY = 0.5 


/i Y = 1.0 



1 


-0.6536170244 
-0.0242156982 


-0.2962340117 


-0.0205712318 



we have used the bound state energies calculated with our third method as an input to 
calculate the wave functions. 



Table 9: Eigenvalues (nonrelativistic) for the case m = 1.0, p = —1.0 and py = 0.02 
calculated by discretizing the integral equation (second method) and making an expansion 
of the eigenfunctions (third method). 



k 


integral eqn. 


exp. in basis 





-0.230696 


-0.2305848598 


1 


-0.044726 


-0.0447072983 


2 


-0.012351 


-0.0123461485 


3 


-0.002981 


-0.0029529333 



Using all three methods, we have calculated the wave functions. In Fig. ^| the wave 
functions in both configuration space and momentum space are given as they are calculated 
by expanding the eigenfunctions in terms of the basis functions. 

We have also calculated the wave function in momentum space by the method of dis- 
cretizing the integral equation (Fig. |Io| ). 

Finally we have calculated the wave functions in coordinate space by directly solving the 
differential equation. We have taken the eigenvalues calculated with our third method as 
input. The wave functions calculated by this method are plotted in Fig. [n]. 

Looking at this figure, we see that the wave functions at the end show some strange 
behaviour, that is they do not converge. The reason that these wave functions do not converge 
is due to the accumulation of errors Note that the wave functions shown in Fig. |ll| are not 
normalized. When the wave functions calculated with both methods are the same their nodes 
and extrema must coincide. In table 10 estimates for hte nodes and extrema are given. 

We see that the zero points of the wave functions lie in the same place, which means that 
we are dealing with the same wave functions. 
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Table 10: Nodes and extrema of the first four wave functions of the nonrelativistic Yukawa 
Hamiltonian calculated with the first method (directly solving the differential equation) and 
the third method (expanding the cigcnfunctions in terms of basis functions). The parameters 
are m — 1.0, p = —1.0 and p,y — 0.02. 





nodes 


extrema 


k 


solving 


diff. eqn. 


exp. in 


basis 


solving 


diff. eqn. 


exp. in 


basis 













1.98 


- 2.01 


2.00 - 


2.01 


1 


3.97 


- 4.04 


4.01 - 


4.02 


1.50 
10.51 


- 1.53 

- 10.61 


1.53 - 
10.56 - 


1.54 
10.58 


2 


3.82 
14.45 


- 3.88 

- 14.59 


3.82 - 
14.50 - 


3.83 
14.51 


1.48 
8.44 
27.28 


- 1.50 

- 8.52 

- 27.50 


1.48 - 

8.49 - 
27.28 - 


1.49 
8.50 
27.33 


3 


3.73 
13.54 
33.45 


- 3.79 

- 13.66 

- 33.72 


3.78 - 
13.64- 
33.46 - 


3.79 

13.65 

33.47 


1.47 
8.16 
22.46 


- 1.50 

- 8.24 

- 22.62 


1.46 - 
8.20 - 
22.41 - 


1.47 
8.21 
22.45 



5 Physical interesting cases 

The physically more interesting situation is the relativistic situation. For the Yukawa po- 
tential we have calculated the bound state energies and wave functions in the relativistic 
situation. In this case we see that the wave functions collapse above a certain strength of the 
potential. This collapse occurs only for certain types of potentials. 

5.1 The relativistic case 

We have done some calculations in the relativistic case for the Yukawa potential. In this 
situation the bound states are calculated with the method of expanding the eigenfunctions 
in terms of basis functions only. We have looked at the same case as in the nonrelativistic 
situation: m = 1.0, p = —1.0 and py = 0.02. The results are given in table |ll|. For 
comparison we have given the nonrelativistic bound states in the same table. The wave 
functions are plotted in Fig. |l2] . 

Table 11: Calculated bound state energies (relativistic and nonrelativistic) for the Yukawa 
potential in the case m = 1.0, p = —1.0 and py = 0.02. 



k 


relativistic 


nonrelativistic 





-0.3107976913 


-0.2305848598 


1 


-0.0583696365 


-0.0447072983 


2 


-0.0161626339 


-0.0123461485 


3 


-0.0042042733 


-0.0029529333 
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5.2 The collapse of wave functions 

It is interesting to look at what happens in the relativistic situation to the bound states and 
the wave functions when the strength (p) of the potential increases. We know that for the 
Coulomb potential there exists a critical value of the strength above which the wave functions 
collapse H We expect to find this collapse for the Yukawa potential as well, since this 
potential is also not bounded from below. 

With a variational calculation we are able to estimate this critical value for the Yukawa 
potential. As the variational wave function we have taken the lowest Olsson wave function 
(which is the same as the lowest Coulomb wave funtion) . These are given in coordinate space 

by 

ip{r) = VVe^ (31) 

and in momentum space by 



= VToTW (32) 

Using these trial functions we can calculate the potential energy exactly: 

r°° p-(iyf 

V(ji) = / drr 2 ^flule'^p VVe^ 



(2ft + fly) 2 ' 
The kinetic energy is given by the integral 



(33) 



7T (p 2 + fl 2 ) 2 V 7T (p 2 + /i 2 ) 



64/x 5 f°° ^ p 2 y/p 2 +m 2 
(p 2 + M 2 ) 



In the limit /i — > oo we can estimate the kinetic energy. Defining = m/fi and p = fj,x and 
substituting this in Eq. (64), we get 



tt y (* 2 + i) 4 

64^ 



,.3 

dx- 



7T 7 (^ + 1) 4 

64/x 
127r' 



(35) 



In the second step we have taken the limit y, — > oo, or /3 — > 0. In this limit, the potential 
energy becomes: 



4p 3 p 
4^ 

For the total Hamiltonian in this limit we find 



^) = TJ=W- ( 36 ) 



^(M) = ^+PM- (37) 
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From this Hamiltonian we can determine the critical value, for which the potential is not 
bounded from below any more. 

Pent = ~ -1.697652726 (38) 

A more accurate estimation is given by ||: p cr n — — 4/-7T = —1.2732395. 

To see whether this collapse really occurs, we have looked at what happens when we 
increase the strength p and keep all other parameters the same. In Fig. [l^ we have plotted 
the wave functions of the ground state in configuration space for m = 1.0, py — 0.02 and 
several values of p. In Fig. [l4] the same is done for the wave functions in the momentum 
space representation. In both figures we see that the wave functions show some irregularities 
at larger p. For p < p cr u the wave function is not a smooth function anymore and ceases to 
be approximated well by our basis. So we get some irregularities due to numerical errors. 

To illustrate that the ground state energy drops very fast, we have plotted it in Fig. |l5| 
for different values of p and several values of py. We see that the curves for different py do 
not coincide, but relative differences are getting smaller for larger p. 

This collapse is typical for potentials that are not bounded from below, like the Coulomb 
and Yukawa potential. To show the difference, we have done the same calculations for the 
Saxon- Woods potential 

Vlr) = , P - w . (39) 

w 1 + e (r-R)/a \ > 

This potential is bounded from below, which means that the bound state energies cannot 
drop to infinity. We have calculated the bound state energies and wave functions in both 
configuration and momentum space for the case 

rn = 1.0, R= 10.0, a = 1.0, p, = 0.35 (40) 

and several values of p. The bound state energies are given in table in Fig. |l| the wave 
functions in configuration space are plotted and in Fig. |l7| they are plotted in the momentum 
space representation. Looking at Figs. |l^ and [ll] we see that the wave functions do not show 
any irregularities. 



Table 12: Bound state energies for the Saxon Woods potential in the relativistic case for 
m = 1.0, R = 10.0, a = 1.0. 



n 


p = -0.1 


p = -0.397975 


p = -1.0 


p = -2.0 





-0.04328906536 


-0.3107855320 


-0.8903754950 


-1.870764375 


1 




-0.1143078804 


-0.6332213879 


-1.569842815 


2 






-0.3176877499 


-1.184906960 


3 






-0.03023195267 


-0.7476220131 



Finally we have plotted the ground state energy versus the strength p for the Saxon 
Woods potential in the same figure with the Yukawa potential (Fig. 15). Comparing these 
two potentials, we notice that the curve of the Saxon Woods potential goes down smoothly 
instead of dropping down, as for the Yukawa potential. The reason for this behaviour is that 
the Saxon Woods potential is bounded from below and does not have a collapse. 
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6 Conclusion 



The easiest way to obtain the bound states is to directly solve the differential equation. This 
method is very easy to implement and calculations do not take up much time. Unfortunately 
it is only uscfull in the nonrelativistic situation. Beside, it is not very accurate. 

In the relativistic situation we use two other methods; discretizing the integral equation 
and making an expansion in terms of basis functions. For both methods the implementation 
is more complicated than in the case of directly solving the differential equation, but it is not 
extremely difficult. Calculations do take up some more time, but this is negligible. And both 
methods are more accurate than the method of directly solving the differential equation. 

A difference between the methods is that the calculations with the method of discretizing 
the integral equation are done in momentum space only. When making an expansion in 
terms of basis functions the kinetic and potential energies are seperately calculated in either 
momentum or configuration space. A consequence of this difference is that the method of 
discretizing the integral equation is not very wellsuited to all potentials, since we may have 
to deal with singularities in the potential matrix elements. 

After checking our codes with the exact solutions of several potentials, we have calculated 
the bound state energies and wave functions for the Yukawa potential in the relativistic 
case. We found a more deeply bound state in the relativistic situation compared to the 
nonrelativistic one. We have also looked at what happens to the wave functions when the 
strength of the potential increases. It turned out that there exists a critical value for the 
strength above which the wave functions collapse. This collapse only occurs for potentials 
that are not bounded from below; e.g. the Coulomb or Yukawa potential. Potentials that 
are bounded from below do not show this collapse, since their ground state energies do not 
drop down to minus infinity. As an example of this kind of potential we have looked at the 
Saxon Woods potential. 
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Figure 1: First five basis function in configuration space (a) and momentum space (b). 




Figure 2: Wave functions in configuration space calculated with (a) 5, (b) 10, (c) 15 or (d) 
20 basis functions for the Yukawa potential in the case m = 1.0, p = —2.0 and /iy = 0.02. 
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Figure 3: Wave functions in momentum space calculated with (a) 5, (b) 10, (c) 15 or (d) 20 
basis functions for the Yukawa potential in the case m = 1.0, p = —2.0 and [iy = 0.02. 




Figure 4: Configuration space wave functions obtained by making an expansion in basis 
functions. Linear potential with m = 1.0, a = 1.0. 
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Figure 5: Calculated wave functions (by making an expansion in basis functions) (a) and exact 
wave functions (b) in configuration space for the Hulthen potential in the case: m = 1.0, 
p = —0.5 and pu = 0.1. 




Figure 6: Calculated wave functions in momentum space for the Hulthen potential for the 
case: m = 1.0, p = —0.5 and pn = 0.1. 
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Figure 7: Wave functions calculated for the Hulthen potential by directly solving the differ- 
ential equation for the case: m = 1.0, p = —0.5 and /xh = 0.1. 
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Figure 8: Bound states of a nonrelativistic qq system in a Yukawa potential for different 
values of fiy. 
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Figure 9: Wave functions in (a) configuration space and (b) momentum space calculated for 
the Yukawa potential by making an expansion of the eigenfunctions in basis functions for the 
case: m = 1.0, p = —1.0 and p\ = 0.02. 




°' 5 x(p=|itan(x)) 1,0 

Figure 10: Wave functions in momentum space calculated for the Yukawa potential by directly 
discretizing the integral equation for the case: m = 1.0, p = —1.0 and py = 0.02. 
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Figure 11: Wave functions in configuration space calculated for the Yukawa potential by 
directly solving the differential equation for the case: m = 1.0, p = —1.0 and py = 0.02. 
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Figure 12: Wave functions for the Yukawa potential by making an expansion in basis functions 
for the case: m = 1.0, p = —1.0 and py — 0.02 in the relativistic situation. 
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Figure 13: Wave function (relativistic) of the ground state in configuration space for 
Yukawa potential with m = 1.0, /iy = 0.02 and several values of p. 
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Figure 14: Wave function (relativistic) of the ground state in momentum space for the Yukawa 
potential with m = 1.0, /iy = 0.02 and several values of p. 
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Figure 16: Wave function (relativistic) of the ground state in configuration space for the 
Saxon Woods potential with m = 1.0, R = 10.0, a = 1.0. 
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Figure 17: Wave function (relativistic) of the ground state in momentum space for the Saxon 
Woods potential with m = 1.0, R = 10.0, a = 1.0. 
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